#include "header.h"

/* --------------------------------------------------------------------- */
void rotationmat1(double x, double y, double z, double R[3][3]){
/* --------------------------------------------------------------------- */
double nx, ny, nz, r;

r = sqrt(x*x + y*y + z*z);
nx= x/r;
ny= y/r;
nz= z/r;

R[0][0] = nx;     R[0][1] = ny;     R[0][2] = nz;
R[1][0] =-ny;     R[1][1] = nx;     R[1][2] = 0.0;
R[2][0] =-nx*nz;  R[2][1] =-ny*nz;  R[2][2] = nx*nx + ny*ny;
}
/* --------------------------------------------------------------------- */
//void rotationmat(double x, double y, double z, double R[3][3]){
void rotationmat(double x, double y, double z, double** R){
/* --------------------------------------------------------------------- */
double phi, sintheta, costheta, sinphi, cosphi, r;

phi = atan2(y, x);
sinphi = sin(phi);
cosphi = cos(phi);

r = sqrt(x*x + y*y + z*z);
sintheta = sin(acos(z/r));
costheta = z/r;

R[0][0] = sintheta*cosphi;  R[0][1] = sintheta*sinphi;  R[0][2] = costheta;
R[1][0] =-sinphi;           R[1][1] = cosphi;           R[1][2] = 0.0;
R[2][0] =-costheta*cosphi;  R[2][1] =-costheta*sinphi;  R[2][2] = sintheta;
}
/* --------------------------------------------------------------------- */
//void rotate(double* var,double R[3][3]){
void rotate(double* var,double** R){
/* --------------------------------------------------------------------- */
double vx,vy,vz;
vx = var[1]; vy = var[2]; vz = var[3];

var[1] = R[0][0]*vx + R[0][1]*vy + R[0][2]*vz;
var[2] = R[1][0]*vx + R[1][1]*vy + R[1][2]*vz;
var[3] = R[2][0]*vx + R[2][1]*vy + R[2][2]*vz;

}
/* --------------------------------------------------------------------- */
//void rotate_back(double* var,double R[3][3]){
void rotate_back(double* var, double** R){
/* --------------------------------------------------------------------- */
double vx,vy,vz;

vx = var[1]; vy = var[2]; vz = var[3];

var[1] = R[0][0]*vx + R[1][0]*vy + R[2][0]*vz;
var[2] = R[0][1]*vx + R[1][1]*vy + R[2][1]*vz;
var[3] = R[0][2]*vx + R[1][2]*vy + R[2][2]*vz;
}
/* --------------------------------------------------------------------- */
void rotate_x(double *x, double *y, double *z, struct node *p){
/* --------------------------------------------------------------------- */
double xx,yy,zz;

xx = *x; yy = *y; zz = *z;

*x    = xx*(p->t1[0]) + yy*(p->t1[1]) + zz*(p->t1[2]);
*y    = xx*(p->t2[0]) + yy*(p->t2[1]) + zz*(p->t2[2]);
*z    = xx*(p->n[0])  + yy*(p->n[1])  + zz*(p->n[2]);
}
/* --------------------------------------------------------------------- */
